Correlations + Mixed models (Igf & animal welfare)

Author

Anja Eggert & Anne-Marie Galow

Published

March 3, 2025

R Libraries

library(tidyverse)     # tidy universe
library(kableExtra)    # html-table
library(patchwork)     # combine plots
library(lmerTest)      # mixed model
library(car)           # ANOVA
library(emmeans)       # post-hoc
library(performance)   # model performance
library(Hmisc)         # Computes correlation and p-value matrix
library(viridis)       # colour scale
set.seed(1989)

Data

Read data

Data processing done in file “igf-biomarker-summary-stats.qmd”.

load("./data/data-processed.RData")

Correlations of parameters

All pairwise correlations

dat.cor.all <- dat.w |> 
  dplyr::select(where(is.numeric))

Calculate correlations and p values:

cor.all <- rcorr(as.matrix(dat.cor.all))

cor.all_matrix  <- round(cor.all$r, 2)
pval.all_matrix <- round(cor.all$P, 3)
cor.all_matrix |> 
  kable(caption = "Person correlation coefficients") |>
  kable_styling(bootstrap_options = c("striped", "hover"), 
                font_size = 12) |> 
  scroll_box()
Person correlation coefficients
insem.age total.piglets prop.males prop.stillborn prop.spreizer prop.lbw mean.bodyweight cort.ser.105dpc cort.ser.8dpp cort.sal.105dpc cort.sal.8dpp bioact.ser.30dpc bioact.ser.105dpc bioact.ser.8dpp igf2.ser.105dpc igf2.ser.8dpp igf1.ser.105dpc igf1.ser.8dpp igfbp2.ser.105dpc igfbp2.ser.8dpp igfbp3.ser.105dpc igfbp3.ser.8dpp proteolysis.105dpc proteolysis.8dpp stc1.105dpc stc1.8dpp calc.30dpc calc.105dpc calc.8dpp igfbp23.ser.105dpc igfbp23.ser.8dpp
insem.age 1.00 0.35 -0.06 0.25 -0.21 -0.20 0.16 -0.19 -0.44 0.02 -0.19 -0.29 -0.20 -0.21 0.07 -0.28 0.53 0.42 0.18 0.32 -0.15 -0.10 0.19 0.13 0.15 0.10 -0.10 0.00 0.01 0.12 0.05
total.piglets 0.35 1.00 0.19 0.33 0.31 0.29 -0.40 0.09 -0.02 0.23 -0.01 0.13 0.05 0.16 -0.14 -0.40 -0.05 -0.04 0.14 0.12 -0.23 -0.11 0.17 -0.30 -0.26 -0.30 0.20 0.32 0.15 0.21 0.05
prop.males -0.06 0.19 1.00 0.09 0.07 0.21 -0.21 0.24 0.49 -0.09 -0.11 -0.03 -0.12 -0.16 0.20 0.18 -0.16 -0.21 0.04 0.18 0.25 0.29 0.08 -0.04 -0.23 -0.22 0.06 0.07 -0.05 -0.16 -0.12
prop.stillborn 0.25 0.33 0.09 1.00 -0.10 -0.11 0.00 -0.09 0.43 -0.06 0.33 -0.22 -0.07 -0.11 0.06 -0.15 -0.19 -0.28 0.08 0.24 0.01 0.09 0.37 0.17 -0.38 -0.40 -0.05 0.09 0.00 -0.11 -0.13
prop.spreizer -0.21 0.31 0.07 -0.10 1.00 0.27 -0.25 0.05 -0.10 0.12 0.05 0.21 0.15 0.27 0.09 -0.24 -0.05 -0.04 0.20 0.09 -0.11 0.05 -0.16 -0.16 -0.44 -0.39 0.30 0.32 -0.15 0.16 0.03
prop.lbw -0.20 0.29 0.21 -0.11 0.27 1.00 -0.88 -0.11 -0.18 0.43 -0.05 0.18 0.10 0.12 0.37 0.28 -0.04 0.28 0.14 -0.12 0.27 -0.03 0.15 -0.42 -0.01 -0.06 -0.07 -0.03 -0.34 0.01 0.03
mean.bodyweight 0.16 -0.40 -0.21 0.00 -0.25 -0.88 1.00 -0.13 0.01 -0.42 -0.04 -0.24 -0.15 -0.13 -0.40 -0.30 -0.05 -0.27 -0.14 0.08 -0.13 0.11 -0.25 0.53 0.00 0.06 0.04 -0.03 0.28 -0.09 -0.10
cort.ser.105dpc -0.19 0.09 0.24 -0.09 0.05 -0.11 -0.13 1.00 0.49 0.31 -0.13 0.24 0.15 0.13 0.26 -0.15 0.23 0.13 0.14 0.35 0.07 0.15 -0.34 0.20 0.06 -0.06 0.23 0.54 0.24 -0.01 0.14
cort.ser.8dpp -0.44 -0.02 0.49 0.43 -0.10 -0.18 0.01 0.49 1.00 -0.22 0.02 -0.40 -0.44 -0.53 0.11 0.09 -0.09 -0.16 0.11 0.40 0.35 0.23 0.02 0.33 -0.08 -0.10 -0.25 0.06 -0.08 -0.24 -0.18
cort.sal.105dpc 0.02 0.23 -0.09 -0.06 0.12 0.43 -0.42 0.31 -0.22 1.00 -0.33 0.51 0.48 0.54 0.40 -0.29 0.33 0.43 0.03 0.13 -0.27 -0.08 0.35 -0.30 0.09 -0.10 0.32 0.24 0.16 0.47 0.32
cort.sal.8dpp -0.19 -0.01 -0.11 0.33 0.05 -0.05 -0.04 -0.13 0.02 -0.33 1.00 -0.35 -0.33 -0.25 -0.15 0.35 -0.03 -0.21 -0.19 -0.48 0.29 0.08 -0.49 0.40 0.04 0.13 0.04 -0.02 -0.05 -0.43 -0.29
bioact.ser.30dpc -0.29 0.13 -0.03 -0.22 0.21 0.18 -0.24 0.24 -0.40 0.51 -0.35 1.00 0.90 0.84 0.00 -0.41 0.21 -0.16 0.11 -0.11 -0.34 -0.27 0.01 -0.34 -0.22 -0.27 0.32 0.34 0.37 0.34 0.37
bioact.ser.105dpc -0.20 0.05 -0.12 -0.07 0.15 0.10 -0.15 0.15 -0.44 0.48 -0.33 0.90 1.00 0.89 0.13 -0.46 0.31 -0.07 0.20 -0.09 -0.29 -0.25 0.07 -0.18 -0.34 -0.38 0.31 0.32 0.26 0.39 0.36
bioact.ser.8dpp -0.21 0.16 -0.16 -0.11 0.27 0.12 -0.13 0.13 -0.53 0.54 -0.25 0.84 0.89 1.00 -0.04 -0.50 0.23 0.02 0.15 -0.11 -0.32 -0.20 0.03 -0.24 -0.25 -0.29 0.50 0.41 0.30 0.43 0.35
igf2.ser.105dpc 0.07 -0.14 0.20 0.06 0.09 0.37 -0.40 0.26 0.11 0.40 -0.15 0.00 0.13 -0.04 1.00 0.13 0.38 0.47 0.50 0.31 0.01 -0.16 0.23 -0.07 -0.13 -0.23 -0.16 -0.14 -0.43 0.44 0.29
igf2.ser.8dpp -0.28 -0.40 0.18 -0.15 -0.24 0.28 -0.30 -0.15 0.09 -0.29 0.35 -0.41 -0.46 -0.50 0.13 1.00 -0.10 0.03 -0.11 -0.15 0.35 0.27 0.00 -0.24 0.53 0.64 -0.21 -0.21 -0.17 -0.41 -0.34
igf1.ser.105dpc 0.53 -0.05 -0.16 -0.19 -0.05 -0.04 -0.05 0.23 -0.09 0.33 -0.03 0.21 0.31 0.23 0.38 -0.10 1.00 0.66 0.09 0.14 -0.36 -0.58 -0.37 -0.29 0.28 0.15 0.06 -0.02 -0.23 0.62 0.87
igf1.ser.8dpp 0.42 -0.04 -0.21 -0.28 -0.04 0.28 -0.27 0.13 -0.16 0.43 -0.21 -0.16 -0.07 0.02 0.47 0.03 0.66 1.00 0.01 -0.04 -0.08 -0.30 -0.23 -0.35 0.39 0.22 0.06 -0.13 -0.35 0.54 0.60
igfbp2.ser.105dpc 0.18 0.14 0.04 0.08 0.20 0.14 -0.14 0.14 0.11 0.03 -0.19 0.11 0.20 0.15 0.50 -0.11 0.09 0.01 1.00 0.68 0.27 0.15 0.10 -0.26 -0.43 -0.46 -0.16 -0.05 -0.27 0.22 0.16
igfbp2.ser.8dpp 0.32 0.12 0.18 0.24 0.09 -0.12 0.08 0.35 0.40 0.13 -0.48 -0.11 -0.09 -0.11 0.31 -0.15 0.14 -0.04 0.68 1.00 0.05 0.22 0.38 -0.23 -0.33 -0.36 -0.21 -0.07 -0.30 0.13 0.13
igfbp3.ser.105dpc -0.15 -0.23 0.25 0.01 -0.11 0.27 -0.13 0.07 0.35 -0.27 0.29 -0.34 -0.29 -0.32 0.01 0.35 -0.36 -0.08 0.27 0.05 1.00 0.72 -0.27 0.24 -0.06 -0.03 -0.13 -0.01 -0.12 -0.62 -0.49
igfbp3.ser.8dpp -0.10 -0.11 0.29 0.09 0.05 -0.03 0.11 0.15 0.23 -0.08 0.08 -0.27 -0.25 -0.20 -0.16 0.27 -0.58 -0.30 0.15 0.22 0.72 1.00 -0.09 0.43 -0.09 -0.02 0.16 0.20 0.05 -0.58 -0.68
proteolysis.105dpc 0.19 0.17 0.08 0.37 -0.16 0.15 -0.25 -0.34 0.02 0.35 -0.49 0.01 0.07 0.03 0.23 0.00 -0.37 -0.23 0.10 0.38 -0.27 -0.09 1.00 -0.30 -0.25 -0.28 -0.31 -0.02 -0.08 0.16 -0.20
proteolysis.8dpp 0.13 -0.30 -0.04 0.17 -0.16 -0.42 0.53 0.20 0.33 -0.30 0.40 -0.34 -0.18 -0.24 -0.07 -0.24 -0.29 -0.35 -0.26 -0.23 0.24 0.43 -0.30 1.00 -0.07 -0.05 0.15 0.08 0.25 -0.54 -0.41
stc1.105dpc 0.15 -0.26 -0.23 -0.38 -0.44 -0.01 0.00 0.06 -0.08 0.09 0.04 -0.22 -0.34 -0.25 -0.13 0.53 0.28 0.39 -0.43 -0.33 -0.06 -0.09 -0.25 -0.07 1.00 0.95 -0.01 -0.16 0.14 -0.12 0.11
stc1.8dpp 0.10 -0.30 -0.22 -0.40 -0.39 -0.06 0.06 -0.06 -0.10 -0.10 0.13 -0.27 -0.38 -0.29 -0.23 0.64 0.15 0.22 -0.46 -0.36 -0.03 -0.02 -0.28 -0.05 0.95 1.00 0.03 -0.13 0.12 -0.23 -0.07
calc.30dpc -0.10 0.20 0.06 -0.05 0.30 -0.07 0.04 0.23 -0.25 0.32 0.04 0.32 0.31 0.50 -0.16 -0.21 0.06 0.06 -0.16 -0.21 -0.13 0.16 -0.31 0.15 -0.01 0.03 1.00 0.67 0.42 0.20 0.04
calc.105dpc 0.00 0.32 0.07 0.09 0.32 -0.03 -0.03 0.54 0.06 0.24 -0.02 0.34 0.32 0.41 -0.14 -0.21 -0.02 -0.13 -0.05 -0.07 -0.01 0.20 -0.02 0.08 -0.16 -0.13 0.67 1.00 0.54 -0.14 -0.11
calc.8dpp 0.01 0.15 -0.05 0.00 -0.15 -0.34 0.28 0.24 -0.08 0.16 -0.05 0.37 0.26 0.30 -0.43 -0.17 -0.23 -0.35 -0.27 -0.30 -0.12 0.05 -0.08 0.25 0.14 0.12 0.42 0.54 1.00 -0.12 -0.15
igfbp23.ser.105dpc 0.12 0.21 -0.16 -0.11 0.16 0.01 -0.09 -0.01 -0.24 0.47 -0.43 0.34 0.39 0.43 0.44 -0.41 0.62 0.54 0.22 0.13 -0.62 -0.58 0.16 -0.54 -0.12 -0.23 0.20 -0.14 -0.12 1.00 0.76
igfbp23.ser.8dpp 0.05 0.05 -0.12 -0.13 0.03 0.03 -0.10 0.14 -0.18 0.32 -0.29 0.37 0.36 0.35 0.29 -0.34 0.87 0.60 0.16 0.13 -0.49 -0.68 -0.20 -0.41 0.11 -0.07 0.04 -0.11 -0.15 0.76 1.00
pval.all_matrix |> 
  kable(caption = "p values of Person correlation coefficients") |>
  kable_styling(bootstrap_options = c("striped", "hover"), 
                font_size = 12) |> 
  scroll_box()
p values of Person correlation coefficients
insem.age total.piglets prop.males prop.stillborn prop.spreizer prop.lbw mean.bodyweight cort.ser.105dpc cort.ser.8dpp cort.sal.105dpc cort.sal.8dpp bioact.ser.30dpc bioact.ser.105dpc bioact.ser.8dpp igf2.ser.105dpc igf2.ser.8dpp igf1.ser.105dpc igf1.ser.8dpp igfbp2.ser.105dpc igfbp2.ser.8dpp igfbp3.ser.105dpc igfbp3.ser.8dpp proteolysis.105dpc proteolysis.8dpp stc1.105dpc stc1.8dpp calc.30dpc calc.105dpc calc.8dpp igfbp23.ser.105dpc igfbp23.ser.8dpp
insem.age NA 0.028 0.719 0.115 0.202 0.207 0.337 0.432 0.049 0.919 0.422 0.075 0.223 0.187 0.754 0.228 0.012 0.050 0.286 0.057 0.379 0.563 0.430 0.578 0.536 0.672 0.530 0.999 0.950 0.468 0.766
total.piglets 0.028 NA 0.230 0.040 0.054 0.071 0.010 0.720 0.920 0.337 0.956 0.442 0.776 0.326 0.545 0.083 0.830 0.846 0.399 0.470 0.178 0.512 0.490 0.193 0.282 0.209 0.221 0.046 0.344 0.214 0.764
prop.males 0.719 0.230 NA 0.589 0.657 0.189 0.200 0.328 0.029 0.728 0.634 0.878 0.467 0.323 0.401 0.457 0.469 0.342 0.834 0.284 0.143 0.090 0.757 0.867 0.336 0.362 0.736 0.662 0.757 0.352 0.489
prop.stillborn 0.115 0.040 0.589 NA 0.555 0.491 0.981 0.720 0.056 0.801 0.153 0.177 0.666 0.494 0.809 0.530 0.392 0.200 0.635 0.151 0.956 0.587 0.116 0.477 0.111 0.091 0.746 0.600 0.985 0.520 0.460
prop.spreizer 0.202 0.054 0.657 0.555 NA 0.091 0.115 0.825 0.687 0.624 0.848 0.209 0.355 0.092 0.705 0.311 0.840 0.873 0.254 0.588 0.538 0.758 0.511 0.496 0.060 0.103 0.068 0.046 0.347 0.360 0.852
prop.lbw 0.207 0.071 0.189 0.491 0.091 NA 0.000 0.664 0.436 0.070 0.842 0.281 0.527 0.473 0.110 0.235 0.857 0.214 0.421 0.477 0.117 0.857 0.550 0.066 0.961 0.820 0.653 0.843 0.030 0.958 0.857
mean.bodyweight 0.337 0.010 0.200 0.981 0.115 0.000 NA 0.582 0.979 0.071 0.858 0.147 0.368 0.421 0.079 0.203 0.829 0.216 0.406 0.644 0.435 0.530 0.306 0.016 0.987 0.795 0.823 0.861 0.084 0.617 0.563
cort.ser.105dpc 0.432 0.720 0.328 0.720 0.825 0.664 0.582 NA 0.032 0.208 0.604 0.323 0.540 0.582 0.300 0.565 0.354 0.590 0.556 0.137 0.785 0.527 0.339 0.565 0.826 0.813 0.353 0.016 0.326 0.974 0.567
cort.ser.8dpp 0.049 0.920 0.029 0.056 0.687 0.436 0.979 0.032 NA 0.366 0.942 0.079 0.055 0.017 0.659 0.712 0.699 0.489 0.649 0.080 0.131 0.326 0.958 0.292 0.746 0.671 0.286 0.812 0.737 0.301 0.446
cort.sal.105dpc 0.919 0.337 0.728 0.801 0.624 0.070 0.071 0.208 0.366 NA 0.172 0.026 0.036 0.017 0.098 0.246 0.168 0.069 0.912 0.601 0.263 0.741 0.325 0.367 0.732 0.692 0.184 0.331 0.524 0.042 0.186
cort.sal.8dpp 0.422 0.956 0.634 0.153 0.848 0.842 0.858 0.604 0.942 0.172 NA 0.132 0.152 0.287 0.552 0.139 0.912 0.376 0.413 0.032 0.221 0.737 0.130 0.200 0.874 0.583 0.882 0.947 0.841 0.058 0.216
bioact.ser.30dpc 0.075 0.442 0.878 0.177 0.209 0.281 0.147 0.323 0.079 0.026 0.132 NA 0.000 0.000 1.000 0.072 0.341 0.485 0.518 0.539 0.043 0.121 0.966 0.143 0.362 0.257 0.050 0.035 0.021 0.046 0.030
bioact.ser.105dpc 0.223 0.776 0.467 0.666 0.355 0.527 0.368 0.540 0.055 0.036 0.152 0.000 NA 0.000 0.577 0.039 0.165 0.761 0.233 0.595 0.085 0.148 0.764 0.436 0.154 0.113 0.055 0.044 0.106 0.018 0.031
bioact.ser.8dpp 0.187 0.326 0.323 0.494 0.092 0.473 0.421 0.582 0.017 0.017 0.287 0.000 0.000 NA 0.882 0.024 0.293 0.920 0.384 0.511 0.057 0.232 0.901 0.303 0.300 0.231 0.001 0.009 0.063 0.009 0.036
igf2.ser.105dpc 0.754 0.545 0.401 0.809 0.705 0.110 0.079 0.300 0.659 0.098 0.552 1.000 0.577 0.882 NA 0.574 0.102 0.037 0.026 0.178 0.966 0.501 0.471 0.820 0.596 0.361 0.499 0.569 0.060 0.054 0.207
igf2.ser.8dpp 0.228 0.083 0.457 0.530 0.311 0.235 0.203 0.565 0.712 0.246 0.139 0.072 0.039 0.024 0.574 NA 0.676 0.895 0.646 0.535 0.127 0.254 1.000 0.430 0.025 0.004 0.381 0.366 0.466 0.074 0.147
igf1.ser.105dpc 0.012 0.830 0.469 0.392 0.840 0.857 0.829 0.354 0.699 0.168 0.912 0.341 0.165 0.293 0.102 0.676 NA 0.001 0.677 0.530 0.105 0.005 0.211 0.310 0.253 0.540 0.806 0.924 0.297 0.002 0.000
igf1.ser.8dpp 0.050 0.846 0.342 0.200 0.873 0.214 0.216 0.590 0.489 0.069 0.376 0.485 0.761 0.920 0.037 0.895 0.001 NA 0.958 0.873 0.723 0.170 0.442 0.225 0.095 0.361 0.793 0.561 0.109 0.009 0.003
igfbp2.ser.105dpc 0.286 0.399 0.834 0.635 0.254 0.421 0.406 0.556 0.649 0.912 0.413 0.518 0.233 0.384 0.026 0.646 0.677 0.958 NA 0.000 0.114 0.367 0.699 0.282 0.066 0.049 0.373 0.785 0.113 0.188 0.341
igfbp2.ser.8dpp 0.057 0.470 0.284 0.151 0.588 0.477 0.644 0.137 0.080 0.601 0.032 0.539 0.595 0.511 0.178 0.535 0.530 0.873 0.000 NA 0.774 0.202 0.124 0.340 0.166 0.131 0.228 0.670 0.077 0.438 0.448
igfbp3.ser.105dpc 0.379 0.178 0.143 0.956 0.538 0.117 0.435 0.785 0.131 0.263 0.221 0.043 0.085 0.057 0.966 0.127 0.105 0.723 0.114 0.774 NA 0.000 0.270 0.325 0.812 0.913 0.456 0.945 0.484 0.000 0.002
igfbp3.ser.8dpp 0.563 0.512 0.090 0.587 0.758 0.857 0.530 0.527 0.326 0.741 0.737 0.121 0.148 0.232 0.501 0.254 0.005 0.170 0.367 0.202 0.000 NA 0.711 0.064 0.715 0.920 0.359 0.251 0.782 0.000 0.000
proteolysis.105dpc 0.430 0.490 0.757 0.116 0.511 0.550 0.306 0.339 0.958 0.325 0.130 0.966 0.764 0.901 0.471 1.000 0.211 0.442 0.699 0.124 0.270 0.711 NA 0.207 0.479 0.442 0.197 0.928 0.741 0.534 0.428
proteolysis.8dpp 0.578 0.193 0.867 0.477 0.496 0.066 0.016 0.565 0.292 0.367 0.200 0.143 0.436 0.303 0.820 0.430 0.310 0.225 0.282 0.340 0.325 0.064 0.207 NA 0.833 0.878 0.524 0.745 0.284 0.018 0.085
stc1.105dpc 0.536 0.282 0.336 0.111 0.060 0.961 0.987 0.826 0.746 0.732 0.874 0.362 0.154 0.300 0.596 0.025 0.253 0.095 0.066 0.166 0.812 0.715 0.479 0.833 NA 0.000 0.960 0.525 0.559 0.614 0.664
stc1.8dpp 0.672 0.209 0.362 0.091 0.103 0.820 0.795 0.813 0.671 0.692 0.583 0.257 0.113 0.231 0.361 0.004 0.540 0.361 0.049 0.131 0.913 0.920 0.442 0.878 0.000 NA 0.909 0.594 0.616 0.341 0.791
calc.30dpc 0.530 0.221 0.736 0.746 0.068 0.653 0.823 0.353 0.286 0.184 0.882 0.050 0.055 0.001 0.499 0.381 0.806 0.793 0.373 0.228 0.456 0.359 0.197 0.524 0.960 0.909 NA 0.000 0.008 0.252 0.808
calc.105dpc 0.999 0.046 0.662 0.600 0.046 0.843 0.861 0.016 0.812 0.331 0.947 0.035 0.044 0.009 0.569 0.366 0.924 0.561 0.785 0.670 0.945 0.251 0.928 0.745 0.525 0.594 0.000 NA 0.000 0.413 0.517
calc.8dpp 0.950 0.344 0.757 0.985 0.347 0.030 0.084 0.326 0.737 0.524 0.841 0.021 0.106 0.063 0.060 0.466 0.297 0.109 0.113 0.077 0.484 0.782 0.741 0.284 0.559 0.616 0.008 0.000 NA 0.472 0.396
igfbp23.ser.105dpc 0.468 0.214 0.352 0.520 0.360 0.958 0.617 0.974 0.301 0.042 0.058 0.046 0.018 0.009 0.054 0.074 0.002 0.009 0.188 0.438 0.000 0.000 0.534 0.018 0.614 0.341 0.252 0.413 0.472 NA 0.000
igfbp23.ser.8dpp 0.766 0.764 0.489 0.460 0.852 0.857 0.563 0.567 0.446 0.186 0.216 0.030 0.031 0.036 0.207 0.147 0.000 0.003 0.341 0.448 0.002 0.000 0.428 0.085 0.664 0.791 0.808 0.517 0.396 0.000 NA

Selected correlations

Select columns and merge data frames:

dat.cor.sel <- dat.w |> 
  dplyr::select(mean.bodyweight, prop.lbw, igf2.ser.105dpc, igf2.ser.8dpp) |> 
  mutate(prop.lbw = prop.lbw/100) |> 
  rename(`Piglet birth weight` = mean.bodyweight,
         `Proportion LBW`      = prop.lbw,
         `IGF2 105dpc`         = igf2.ser.105dpc,
         `IGF2 8dpp`           = igf2.ser.8dpp)

Calculate correlations and p values:

cor.sel <- rcorr(as.matrix(dat.cor.sel))

cor.sel_matrix  <- as_tibble(cor.sel$r, rownames = "Var1")
pval.sel_matrix <- as_tibble(cor.sel$P, rownames = "Var1")
cor.sel_matrix |> 
  kable(caption = "Person correlation coefficients") |>
  kable_styling(bootstrap_options = c("striped", "hover"), 
                font_size = 12) |> 
  scroll_box()
Person correlation coefficients
Var1 Piglet birth weight Proportion LBW IGF2 105dpc IGF2 8dpp
Piglet birth weight 1.0000000 -0.8823905 -0.4013240 -0.2974563
Proportion LBW -0.8823905 1.0000000 0.3686377 0.2780085
IGF2 105dpc -0.4013240 0.3686377 1.0000000 0.1336761
IGF2 8dpp -0.2974563 0.2780085 0.1336761 1.0000000
pval.sel_matrix |> 
  kable(caption = "p values of Person correlation coefficients") |>
  kable_styling(bootstrap_options = c("striped", "hover"), 
                font_size = 12) |> 
  scroll_box()
p values of Person correlation coefficients
Var1 Piglet birth weight Proportion LBW IGF2 105dpc IGF2 8dpp
Piglet birth weight NA 0.0000000 0.0794637 0.2027779
Proportion LBW 0.0000000 NA 0.1097394 0.2353017
IGF2 105dpc 0.0794637 0.1097394 NA 0.5742138
IGF2 8dpp 0.2027779 0.2353017 0.5742138 NA

Heatmap:

# Convert matrices into long format for ggplot
cor_long <- pivot_longer(cor.sel_matrix,
                         cols = !Var1,
                         names_to = "Var2",
                         values_to = "cor.val")

p_long   <- pivot_longer(pval.sel_matrix,
                         cols = !Var1,
                         names_to = "Var2",
                         values_to = "p.val")

# Merge correlation values and p-values
heatmap_data <- cor_long |> 
  left_join(p_long, by = c("Var1", "Var2")) |> 
  # this makes it to characters
  mutate(cor.val = sprintf("%.2f", cor.val),
         p.val = sprintf("%.3f", p.val)) |> 
  mutate(cor.val = case_when(cor.val == "1.00" ~ NA_character_, 
                             TRUE ~ cor.val)) |> 
  mutate(p.val = case_when(p.val == "0.000" ~ "<0.001",
                           is.na(p.val) ~ NA_character_, 
                           TRUE ~ paste0("=",p.val))) |>
  mutate(label = case_when(is.na(cor.val) ~ NA_character_, 
                           TRUE ~ paste0(cor.val, "\n(p", p.val, ")")))
# Create the heatmap
plot.corr <- heatmap_data |> 
  ggplot(aes(Var1, Var2, fill = as.numeric(cor.val))) +
  geom_tile(color = "white") +
  scale_fill_viridis(limits = c(-1, 1), guide="none") +
  # Overlay correlation & p-values
  geom_text(aes(label = label), color = "white", size = 4) + 
  my_theme_cor +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "",
       x = "",
       y = "")
plot.corr

Statistical model: general design

To correctly account for the nested data structure in the analysis using a linear mixed model with the lmerTest package in R, we need to account for the hierarchical structure in the data:

  • Husbandry Level (conventional/ecological): Each sow is housed in either of the two facilities.
  • Sow Level: Each sow gave birth 1, 2 or 3 times, so these events are nested within sows.
  • Time Point Level: Measurements were taken before and after each birth, introducing repeated measures within each pregnancy.
  • control for age of the sows

This type of nested structure is best represented by a model with random effects at the different levels to account for the dependencies within each level. The random effects allow for variation within subjects and

contr = lmerControl(optimizer   = "bobyqa",
                    optCtrl     = list(maxfun = 10000000),
                    calc.derivs = FALSE)

Cortisol serum

Model

hist(dat.l |> 
       dplyr::filter(parameter == "cort.ser") |> 
       droplevels() |>
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
         dplyr::filter(parameter == "cort.ser") |> 
         droplevels() |>
         drop_na(value) |> 
         pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.cort1 <- lmerTest::lmer(log(value) ~ 
                              husbandry +
                              time +
                              # interaction term
                              husbandry : time +
                              # random intercept for sows and for each litter within sow
                              (1 | sow/litter.no),
                            data    = dat.l |> 
                              dplyr::filter(parameter == "cort.ser") |> 
                              droplevels() |>
                              drop_na(value),
                            REML    = TRUE,
                            control = contr)
summary(mod.cort1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "cort.ser")),  
    value)
Control: contr

REML criterion at convergence: 22.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.45942 -0.49449 -0.08495  0.57811  1.38490 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.02937  0.1714  
 sow           (Intercept) 0.01239  0.1113  
 Residual                  0.05441  0.2333  
Number of obs: 39, groups:  litter.no:sow, 20; sow, 14

Fixed effects:
                             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                   3.40152    0.10154 14.82783  33.498 2.18e-15 ***
husbandryecological          -0.03632    0.14635 16.13178  -0.248    0.807    
time8dpp                     -0.10436    0.10432 17.29836  -1.000    0.331    
husbandryecological:time8dpp  0.11841    0.15039 17.69062   0.787    0.442    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.694              
time8dpp    -0.514  0.356       
hsbndrycl:8  0.356 -0.532 -0.694
round(drop1(mod.cort1, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
husbandry:time  0.034   0.034     1 17.691    0.62  0.442

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.cort1,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                Chisq Df Pr(>Chisq)
husbandry      0.0406  1     0.8404
time           0.3976  1     0.5283
husbandry:time 0.6199  1     0.4311
car::Anova(mod.cort1,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                    F Df Df.res Pr(>F)
husbandry      0.0388  1 10.877 0.8474
time           0.4022  1 17.439 0.5342
husbandry:time 0.6151  1 17.473 0.4434

Model diagnostics

performance::check_model(mod.cort1)

Emmeans & Effect sizes

Most popular effect-size measure is probably Cohen’s d, which is defined as the observed difference, divided by the population SD; and obviously Cohen effect sizes are close cousins of pairwise differences. They are available via the eff_size() function of the emmeans package.

Husbandry

Emmeans:

emm <- emmeans(mod.cort1,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "cort.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response   SE   df lower.CL upper.CL
 conventional     28.5 2.48 8.61     23.4     34.7
 ecological       29.1 2.57 9.27     23.9     35.5

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological 0.977 0.121 8.93    1  -0.185  0.8577

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.cort1),      # Residual standard deviation from model
         edf = df.residual(mod.cort1))  # Residual degrees of freedom
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)  -0.0981 0.532 8.93     -1.3     1.11

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.2333 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.cort1,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "cort.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response   SE   df lower.CL upper.CL
 105dpc     29.5 2.16 16.1     25.2     34.4
 8dpp       28.2 2.02 14.8     24.2     32.8

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio     SE   df null t.ratio p.value
 105dpc / 8dpp  1.05 0.0787 17.7    1   0.600  0.5558

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.cort1),
         edf = df.residual(mod.cort1))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE   df lower.CL upper.CL
 (105dpc - 8dpp)    0.194 0.323 17.7   -0.486    0.874

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.2333 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.cort1,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "cort.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response   SE   df lower.CL upper.CL
 conventional     30.0 3.05 14.8     24.2     37.3
 ecological       28.9 3.05 17.4     23.2     36.1

time = 8dpp:
 husbandry    response   SE   df lower.CL upper.CL
 conventional     27.0 2.75 14.8     21.8     33.6
 ecological       29.3 2.98 14.8     23.6     36.5

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological 1.037 0.152 16.1    1   0.248  0.8071

time = 8dpp:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological 0.921 0.132 14.8    1  -0.572  0.5762

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.cort1),
         edf = df.residual(mod.cort1))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)    0.156 0.628 16.1    -1.17    1.485

time = 8dpp:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)   -0.352 0.617 14.8    -1.67    0.965

sigma used for effect sizes: 0.2333 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.cort1,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "cort.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response   SE   df lower.CL upper.CL
 105dpc     30.0 3.05 14.8     24.2     37.3
 8dpp       27.0 2.75 14.8     21.8     33.6

husbandry = ecological:
 time   response   SE   df lower.CL upper.CL
 105dpc     28.9 3.05 17.4     23.2     36.1
 8dpp       29.3 2.98 14.8     23.6     36.5

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio    SE   df null t.ratio p.value
 105dpc / 8dpp 1.110 0.116 17.3    1   1.000  0.3309

husbandry = ecological:
 contrast      ratio    SE   df null t.ratio p.value
 105dpc / 8dpp 0.986 0.107 18.1    1  -0.130  0.8982

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.cort1),
         edf = df.residual(mod.cort1))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE   df lower.CL upper.CL
 (105dpc - 8dpp)   0.4474 0.451 17.3   -0.502    1.397

husbandry = ecological:
 contrast        estimate    SE   df lower.CL upper.CL
 (105dpc - 8dpp)  -0.0602 0.464 18.1   -1.036    0.915

sigma used for effect sizes: 0.2333 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.cort1 <- dat.l  |> 
  dplyr::filter(parameter == "cort.ser") |> 
  droplevels() |>
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 60), 
                     breaks = seq(0, 60, 20),
                     minor_breaks = seq(0, 60, by = 2) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "Cortisol (serum) [ng/ml]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.cort1

Cortisol saliva

Model

hist(dat.l |> 
       dplyr::filter(parameter == "cort.sal") |> 
       droplevels() |>
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
         dplyr::filter(parameter == "cort.sal") |> 
         droplevels() |>
         drop_na(value) |> 
         pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.cort2 <- lmerTest::lmer(log(value) ~ 
                              husbandry +
                              time +
                              # interaction term
                              husbandry : time +
                              # random intercept for sows and for each litter within sow
                              (1 | sow/litter.no),
                            data    = dat.l |> 
                              dplyr::filter(parameter == "cort.sal") |> 
                              droplevels() |>
                              drop_na(value),
                            REML    = TRUE,
                            control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.cort2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "cort.sal")),  
    value)
Control: contr

REML criterion at convergence: 101.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.29670 -0.33630  0.06505  0.77409  1.36353 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.000    0.0000  
 sow           (Intercept) 0.000    0.0000  
 Residual                  0.827    0.9094  
Number of obs: 39, groups:  litter.no:sow, 20; sow, 14

Fixed effects:
                             Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)                    2.0473     0.3031 35.0000   6.754 7.96e-08 ***
husbandryecological           -0.4572     0.4179 35.0000  -1.094   0.2813    
time8dpp                      -0.7537     0.4179 35.0000  -1.804   0.0799 .  
husbandryecological:time8dpp   0.9303     0.5831 35.0000   1.595   0.1196    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.725              
time8dpp    -0.725  0.526       
hsbndrycl:8  0.520 -0.717 -0.717
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.cort2, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time  2.105   2.105     1    35   2.546   0.12

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.cort2,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                Chisq Df Pr(>Chisq)
husbandry      0.0049  1     0.9439
time           0.8969  1     0.3436
husbandry:time 2.5455  1     0.1106
car::Anova(mod.cort2,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                    F Df  Df.res Pr(>F)
husbandry      0.0058  1  9.7306 0.9409
time           0.8843  1 17.6972 0.3597
husbandry:time 2.5370  1 17.7258 0.1289

Model diagnostics

performance::check_model(mod.cort2)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.cort2,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "cort.sal") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response   SE df lower.CL upper.CL
 conventional     5.31 1.11 35     3.48     8.12
 ecological       5.36 1.09 35     3.55     8.09

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio    SE df null t.ratio p.value
 conventional / ecological 0.992 0.289 35    1  -0.027  0.9785

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.cort2),
         edf = df.residual(mod.cort2))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological) -0.00871 0.321 35    -0.66    0.642

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.9094 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.cort2,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "cort.sal") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response    SE df lower.CL upper.CL
 105dpc     6.16 1.290 35     4.03     9.42
 8dpp       4.62 0.939 35     3.06     6.98

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp  1.33 0.389 35    1   0.990  0.3291

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.cort2),
         edf = df.residual(mod.cort2))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    0.317 0.323 35   -0.338    0.973

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.9094 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.cort2,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "cort.sal") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response   SE df lower.CL upper.CL
 conventional     7.75 2.35 35     4.19    14.34
 ecological       4.90 1.41 35     2.74     8.79

time = 8dpp:
 husbandry    response   SE df lower.CL upper.CL
 conventional     3.65 1.05 35     2.03     6.54
 ecological       5.85 1.68 35     3.26    10.49

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio    SE df null t.ratio p.value
 conventional / ecological 1.580 0.660 35    1   1.094  0.2813

time = 8dpp:
 contrast                  ratio    SE df null t.ratio p.value
 conventional / ecological 0.623 0.253 35    1  -1.163  0.2526

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.cort2),
         edf = df.residual(mod.cort2))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)    0.503 0.464 35   -0.439    1.444

time = 8dpp:
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)   -0.520 0.452 35   -1.438    0.397

sigma used for effect sizes: 0.9094 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.cort2,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "cort.sal") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response   SE df lower.CL upper.CL
 105dpc     7.75 2.35 35     4.19    14.34
 8dpp       3.65 1.05 35     2.03     6.54

husbandry = ecological:
 time   response   SE df lower.CL upper.CL
 105dpc     4.90 1.41 35     2.74     8.79
 8dpp       5.85 1.68 35     3.26    10.49

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp 2.125 0.888 35    1   1.804  0.0799

husbandry = ecological:
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp 0.838 0.341 35    1  -0.434  0.6668

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.cort2),
         edf = df.residual(mod.cort2))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    0.829 0.471 35   -0.127    1.785

husbandry = ecological:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)   -0.194 0.448 35   -1.103    0.715

sigma used for effect sizes: 0.9094 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.cort2 <- dat.l  |> 
  dplyr::filter(parameter == "cort.sal") |> 
  droplevels() |>
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 30), 
                     breaks = seq(0, 30, 10),
                     minor_breaks = seq(0, 30, by = 2) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "Cortisol (saliva) [ng/mg protein]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.cort2

IGF bioactivity

Model

hist(dat.l |> 
       dplyr::filter(parameter == "bioact.ser") |> 
       droplevels() |>
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
         dplyr::filter(parameter == "bioact.ser") |> 
         droplevels() |>
         drop_na(value) |> 
         pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.bioact <- lmerTest::lmer(log(value) ~ 
                               husbandry +
                               time +
                               # interaction term
                               husbandry : time +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "bioact.ser") |> 
                               drop_na(value) |> 
                               dplyr::filter(time != "30dpc") |> 
                               droplevels(),
                            REML    = TRUE,
                            control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.bioact)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: droplevels(dplyr::filter(drop_na(dplyr::filter(dat.l, parameter ==  
    "bioact.ser"), value), time != "30dpc"))
Control: contr

REML criterion at convergence: 143.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.29365 -0.36326  0.07091  0.43812  1.94984 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.3260   0.5709  
 sow           (Intercept) 0.0000   0.0000  
 Residual                  0.1372   0.3703  
Number of obs: 80, groups:  litter.no:sow, 40; sow, 15

Fixed effects:
                             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                   5.11600    0.15217 50.82288  33.620  < 2e-16 ***
husbandryecological          -0.76647    0.21520 50.82288  -3.562 0.000811 ***
time8dpp                     -0.05793    0.11711 38.00000  -0.495 0.623690    
husbandryecological:time8dpp  0.20329    0.16562 38.00000   1.227 0.227228    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.707              
time8dpp    -0.385  0.272       
hsbndrycl:8  0.272 -0.385 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.bioact, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time  0.207   0.207     1    38   1.506  0.227

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.bioact,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                 Chisq Df Pr(>Chisq)    
husbandry      11.2029  1  0.0008167 ***
time            0.2786  1  0.5976150    
husbandry:time  1.5065  1  0.2196781    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.bioact,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                     F Df Df.res  Pr(>F)   
husbandry      11.0519  1 11.923 0.00611 **
time            0.2786  1 38.000 0.60068   
husbandry:time  1.5065  1 38.000 0.22723   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.bioact)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.bioact,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "bioact.ser") |> 
          drop_na(value) |> 
          dplyr::filter(time != "30dpc") |> 
          droplevels(), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response   SE df lower.CL upper.CL
 conventional    161.9 22.7 38    121.8      215
 ecological       83.3 11.7 38     62.7      111

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio    SE df null t.ratio p.value
 conventional / ecological  1.94 0.386 38    1   3.347  0.0019

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.bioact),
         edf = df.residual(mod.bioact))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)      1.8 0.557 38    0.669     2.92

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.3703 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.cort2,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "bioact.ser") |> 
          drop_na(value) |> 
          dplyr::filter(time != "30dpc") |> 
          droplevels(), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response    SE df lower.CL upper.CL
 105dpc     6.16 1.290 35     4.03     9.42
 8dpp       4.62 0.939 35     3.06     6.98

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp  1.33 0.389 35    1   0.990  0.3291

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.bioact),
         edf = df.residual(mod.bioact))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate   SE df lower.CL upper.CL
 (105dpc - 8dpp)    0.779 0.79 35   -0.824     2.38

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.3703 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.bioact,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "bioact.ser") |> 
          drop_na(value) |> 
          dplyr::filter(time != "30dpc") |> 
          droplevels(), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response   SE   df lower.CL upper.CL
 conventional    166.7 25.4 50.8    122.8      226
 ecological       77.4 11.8 50.8     57.1      105

time = 8dpp:
 husbandry    response   SE   df lower.CL upper.CL
 conventional    157.3 23.9 50.8    115.9      213
 ecological       89.6 13.6 50.8     66.0      122

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological  2.15 0.463 50.8    1   3.562  0.0008

time = 8dpp:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological  1.76 0.378 50.8    1   2.617  0.0117

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.bioact),
         edf = df.residual(mod.bioact))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)     2.07 0.606 50.8    0.853     3.29

time = 8dpp:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)     1.52 0.595 50.8    0.327     2.71

sigma used for effect sizes: 0.3703 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.bioact,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "bioact.ser") |> 
          drop_na(value) |> 
          dplyr::filter(time != "30dpc") |> 
          droplevels(), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response   SE   df lower.CL upper.CL
 105dpc    166.7 25.4 50.8    122.8      226
 8dpp      157.3 23.9 50.8    115.9      213

husbandry = ecological:
 time   response   SE   df lower.CL upper.CL
 105dpc     77.4 11.8 50.8     57.1      105
 8dpp       89.6 13.6 50.8     66.0      122

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp 1.060 0.124 38    1   0.495  0.6237

husbandry = ecological:
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp 0.865 0.101 38    1  -1.241  0.2222

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.bioact),
         edf = df.residual(mod.bioact))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    0.156 0.316 38   -0.484    0.797

husbandry = ecological:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)   -0.392 0.318 38   -1.036    0.251

sigma used for effect sizes: 0.3703 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.bioact <- dat.l  |> 
  dplyr::filter(parameter == "bioact.ser") |> 
  drop_na(value) |> 
  dplyr::filter(time != "30dpc") |> 
  droplevels() |>
  # plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 600), 
                     breaks = seq(0, 600, 200),
                     minor_breaks = seq(0, 600, by = 50) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "IGF bioactivity [ng/ml]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.bioact

IGF1 (serum)

Model

hist(dat.l |> 
       dplyr::filter(parameter == "igf1.ser") |> 
       droplevels() |> 
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
       dplyr::filter(parameter == "igf1.ser") |> 
       droplevels() |> 
       drop_na(value) |> 
       pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.igf1 <- lmerTest::lmer(log(value) ~ 
                              husbandry +
                              time +
                              # interaction term
                              husbandry : time +
                              # random intercept for sows and for each litter within sow
                              (1 | sow/litter.no),
                            data    = dat.l |> 
                              dplyr::filter(parameter == "igf1.ser") |> 
                              droplevels() |>
                              drop_na(value),
                            REML    = TRUE,
                            control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.igf1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igf1.ser")),  
    value)
Control: contr

REML criterion at convergence: 77

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.90876 -0.52177 -0.03551  0.52631  1.57329 

Random effects:
 Groups        Name        Variance  Std.Dev. 
 litter.no:sow (Intercept) 2.225e-01 4.717e-01
 sow           (Intercept) 6.468e-17 8.042e-09
 Residual                  1.642e-01 4.052e-01
Number of obs: 44, groups:  litter.no:sow, 22; sow, 15

Fixed effects:
                             Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)                    4.8938     0.1966 30.0492  24.887   <2e-16 ***
husbandryecological           -0.5121     0.2663 30.0492  -1.923   0.0640 .  
time8dpp                       0.5069     0.1812 20.0000   2.797   0.0111 *  
husbandryecological:time8dpp   0.4368     0.2453 20.0000   1.780   0.0902 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.739              
time8dpp    -0.461  0.340       
hsbndrycl:8  0.340 -0.461 -0.739
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.igf1, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
husbandry:time   0.52    0.52     1    20    3.17   0.09 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.igf1,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                 Chisq Df Pr(>Chisq)    
husbandry       1.5445  1    0.21395    
time           37.2036  1  1.064e-09 ***
husbandry:time  3.1696  1    0.07502 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.igf1,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                     F Df Df.res    Pr(>F)    
husbandry       1.4352  1 11.177   0.25570    
time           37.2036  1 20.000 5.826e-06 ***
husbandry:time  3.1696  1 20.000   0.09022 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.igf1)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.igf1,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igf1.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response   SE df lower.CL upper.CL
 conventional      172 30.0 20      119      247
 ecological        128 20.4 20       92      179

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio    SE df null t.ratio p.value
 conventional / ecological  1.34 0.317 20    1   1.243  0.2283

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igf1),
         edf = df.residual(mod.igf1))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)    0.725 0.589 20   -0.504     1.95

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.4052 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.igf1,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igf1.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response   SE   df lower.CL upper.CL
 105dpc      103 13.8 30.1     78.7      136
 8dpp        213 28.4 30.1    162.6      280

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.484 0.0594 20    1  -5.912  <.0001

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igf1),
         edf = df.residual(mod.igf1))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    -1.79 0.367 20    -2.56    -1.02

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.4052 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.igf1,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igf1.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response   SE   df lower.CL upper.CL
 conventional      133 26.2 30.1     89.3      199
 ecological         80 14.4 30.1     55.4      115

time = 8dpp:
 husbandry    response   SE   df lower.CL upper.CL
 conventional      222 43.6 30.1    148.3      331
 ecological        205 36.9 30.1    142.4      296

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological  1.67 0.444 30.1    1   1.923  0.0640

time = 8dpp:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological  1.08 0.287 30.1    1   0.283  0.7793

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igf1),
         edf = df.residual(mod.igf1))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)    1.264 0.673 30.1   -0.111     2.64

time = 8dpp:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)    0.186 0.657 30.1   -1.157     1.53

sigma used for effect sizes: 0.4052 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.igf1,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igf1.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response   SE   df lower.CL upper.CL
 105dpc      133 26.2 30.1     89.3      199
 8dpp        222 43.6 30.1    148.3      331

husbandry = ecological:
 time   response   SE   df lower.CL upper.CL
 105dpc       80 14.4 30.1     55.4      115
 8dpp        205 36.9 30.1    142.4      296

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.602 0.1090 20    1  -2.797  0.0111

husbandry = ecological:
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.389 0.0644 20    1  -5.705  <.0001

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igf1),
         edf = df.residual(mod.igf1))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate   SE df lower.CL upper.CL
 (105dpc - 8dpp)    -1.25 0.47 20    -2.23    -0.27

husbandry = ecological:
 contrast        estimate   SE df lower.CL upper.CL
 (105dpc - 8dpp)    -2.33 0.49 20    -3.35    -1.31

sigma used for effect sizes: 0.4052 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.igf1 <- dat.l  |> 
  dplyr::filter(parameter == "igf1.ser") |> 
  droplevels() |>
  drop_na(value) |> 
  # plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 500), 
                     breaks = seq(0, 500, 100),
                     minor_breaks = seq(0, 500, by = 20) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "IGF1 (serum) [ng/ml]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.igf1

IGF2 (serum)

Model

hist(dat.l |> 
       dplyr::filter(parameter == "igf2.ser") |> 
       droplevels() |> 
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
       dplyr::filter(parameter == "igf2.ser") |> 
       droplevels() |>
       drop_na(value) |> 
       pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.igf2 <- lmerTest::lmer(log(value) ~ 
                              husbandry +
                              time +
                              # interaction term
                              husbandry : time +
                              # random intercept for sows and for each litter within sow
                              (1 | sow/litter.no),
                            data    = dat.l |> 
                              dplyr::filter(parameter == "igf2.ser") |> 
                              droplevels() |> 
                              drop_na(value),
                            REML    = TRUE,
                            control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.igf2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igf2.ser")),  
    value)
Control: contr

REML criterion at convergence: 15.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3607 -0.5221  0.0581  0.4752  1.6754 

Random effects:
 Groups        Name        Variance  Std.Dev. 
 litter.no:sow (Intercept) 1.265e-02 1.125e-01
 sow           (Intercept) 5.401e-17 7.349e-09
 Residual                  5.793e-02 2.407e-01
Number of obs: 40, groups:  litter.no:sow, 20; sow, 14

Fixed effects:
                             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                   4.26837    0.08401 34.87974  50.808  < 2e-16 ***
husbandryecological          -0.33300    0.11881 34.87974  -2.803 0.008215 ** 
time8dpp                      0.42980    0.10764 18.00000   3.993 0.000853 ***
husbandryecological:time8dpp  0.43848    0.15222 18.00000   2.881 0.009955 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.707              
time8dpp    -0.641  0.453       
hsbndrycl:8  0.453 -0.641 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.igf2, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)   
husbandry:time  0.481   0.481     1    18   8.298   0.01 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.igf2,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                 Chisq Df Pr(>Chisq)    
husbandry       1.5551  1    0.21239    
time           72.7197  1    < 2e-16 ***
husbandry:time  8.2975  1    0.00397 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.igf2,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                     F Df Df.res    Pr(>F)    
husbandry       1.4218  1 10.105  0.260359    
time           72.7197  1 18.000 9.755e-08 ***
husbandry:time  8.2975  1 18.000  0.009955 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.igf2)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.igf2,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igf2.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response   SE df lower.CL upper.CL
 conventional     88.5 5.71 18     77.3    101.4
 ecological       79.0 5.10 18     69.0     90.5

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio    SE df null t.ratio p.value
 conventional / ecological  1.12 0.102 18    1   1.247  0.2284

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igf2),
         edf = df.residual(mod.igf2))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)    0.473 0.383 18   -0.333     1.28

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.2407 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.igf2,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igf2.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response   SE   df lower.CL upper.CL
 105dpc     60.5 3.59 34.9     53.6     68.2
 8dpp      115.7 6.87 34.9    102.5    130.5

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.523 0.0398 18    1  -8.528  <.0001

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igf2),
         edf = df.residual(mod.igf2))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)     -2.7 0.458 18    -3.66    -1.73

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.2407 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.igf2,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igf2.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response    SE   df lower.CL upper.CL
 conventional     71.4  6.00 34.9     60.2     84.7
 ecological       51.2  4.30 34.9     43.2     60.7

time = 8dpp:
 husbandry    response    SE   df lower.CL upper.CL
 conventional    109.7  9.22 34.9     92.5    130.2
 ecological      122.0 10.20 34.9    102.8    144.6

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological   1.4 0.166 34.9    1   2.803  0.0082

time = 8dpp:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological   0.9 0.107 34.9    1  -0.888  0.3807

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igf2),
         edf = df.residual(mod.igf2))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)    1.384 0.522 34.9    0.323     2.44

time = 8dpp:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)   -0.438 0.497 34.9   -1.446     0.57

sigma used for effect sizes: 0.2407 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.igf2,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igf2.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response    SE   df lower.CL upper.CL
 105dpc     71.4  6.00 34.9     60.2     84.7
 8dpp      109.7  9.22 34.9     92.5    130.2

husbandry = ecological:
 time   response    SE   df lower.CL upper.CL
 105dpc     51.2  4.30 34.9     43.2     60.7
 8dpp      122.0 10.20 34.9    102.8    144.6

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.651 0.0700 18    1  -3.993  0.0009

husbandry = ecological:
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.420 0.0452 18    1  -8.067  <.0001

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igf2),
         edf = df.residual(mod.igf2))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    -1.79 0.498 18    -2.83   -0.739

husbandry = ecological:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    -3.61 0.630 18    -4.93   -2.284

sigma used for effect sizes: 0.2407 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.igf2 <- dat.l  |> 
  dplyr::filter(parameter == "igf2.ser") |> 
  droplevels() |> 
  drop_na(value) |> 
  # plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 200), 
                     breaks = seq(0, 200, 50),
                     minor_breaks = seq(0, 200, by = 10) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "IGF2 (serum) [ng/ml]",
       fill = "Husbandry") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5))) +
  my_theme +
  theme(legend.position = c(0.25, 0.85),
        legend.box = "horizontal")
Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
plot.igf2

Combined plots: Figure 3

# Combine plots with a designated area for the legend
combined <- (plot.cort1 + 
             plot.cort2 + 
             plot.bioact + 
             plot.igf1) + 
  plot_layout(guides = "collect") +
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(face = "bold", size = 20),
        legend.position = "bottom",
        legend.direction = "horizontal")
combined

png("./plots/figure3.png",
     width = 300, height = 300, units = "mm",
     pointsize = 10, res = 600)

combined

dev.off()
png 
  2 

Combined plots: Figure 4

# Combine plots with a designated area for the legend
combined <- (plot.igf2 +
             plot.corr) + 
  #plot_layout(guides = "collect") +
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(face = "bold", size = 20))
combined

png("./plots/figure4.png",
     width = 300, height = 200, units = "mm",
     pointsize = 10, res = 600)

combined
Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_text()`).
dev.off()
png 
  2 

Litter numbers

plot.bio.li <- dat.w  |> 
  # plot
  mutate(jit = jitter(as.numeric(litter.no), 0.3)) |>  
  ggplot(aes(y   = bioact.ser.105dpc)) +
  geom_boxplot(aes(x   = litter.no, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = litter.no,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 600), 
                     breaks = seq(0, 600, 200),
                     minor_breaks = seq(0, 600, by = 50) ) +
  scale_x_discrete(labels= c("1st", "2nd", "3rd")) +
  labs(x = "Litter number",
       y = "IGF bioactivity 105dpc [ng/ml]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "bottom") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.lbw.li <- dat.w  |> 
  # plot
  mutate(jit = jitter(as.numeric(litter.no), 0.3)) |>  
  ggplot(aes(y   = prop.lbw)) +
  geom_boxplot(aes(x   = litter.no, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = litter.no,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 100), 
                     breaks = seq(0, 100, 20),
                     minor_breaks = seq(0, 100, by = 5) ) +
  scale_x_discrete(labels= c("1st", "2nd", "3rd")) +
  labs(x = "Litter number",
       y = "Proportion LBW [%]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "bottom") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))

Combined plots: Figure 5

# Combine plots with a designated area for the legend
combined <- (plot.bio.li + 
             plot.lbw.li) + 
  plot_layout(guides = "collect") +
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(face = "bold", size = 20),
        legend.position = "bottom")
combined

png("./plots/figure5.png",
     width = 300, height = 200, units = "mm",
     pointsize = 10, res = 600)

combined

dev.off()
png 
  2 

IGFBP2 (serum)

Model

hist(dat.l |> 
       dplyr::filter(parameter == "igfbp2.ser") |> 
       droplevels() |> 
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
         dplyr::filter(parameter == "igfbp2.ser") |> 
         droplevels() |>
         drop_na(value) |> 
         pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.igfbp2 <- lmerTest::lmer(log(value) ~ 
                               husbandry +
                               time +
                               # interaction term
                               husbandry : time +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "igfbp2.ser") |> 
                               droplevels() |> 
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
summary(mod.igfbp2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp2.ser")),  
    value)
Control: contr

REML criterion at convergence: 63.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.80571 -0.56601  0.02105  0.51087  1.84289 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.07463  0.2732  
 sow           (Intercept) 0.04330  0.2081  
 Residual                  0.05864  0.2422  
Number of obs: 72, groups:  litter.no:sow, 36; sow, 15

Fixed effects:
                             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                   6.51548    0.11794 14.52980  55.243   <2e-16 ***
husbandryecological          -0.46275    0.16544 15.40369  -2.797   0.0133 *  
time8dpp                     -0.03649    0.08072 34.00000  -0.452   0.6541    
husbandryecological:time8dpp  0.13669    0.11416 34.00000   1.197   0.2394    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.713              
time8dpp    -0.342  0.244       
hsbndrycl:8  0.242 -0.345 -0.707
round(drop1(mod.igfbp2, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time  0.084   0.084     1    34   1.434  0.239

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.igfbp2,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                Chisq Df Pr(>Chisq)  
husbandry      6.4509  1    0.01109 *
time           0.3114  1    0.57680  
husbandry:time 1.4337  1    0.23116  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.igfbp2,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                    F Df Df.res  Pr(>F)  
husbandry      6.3339  1  12.02 0.02704 *
time           0.3114  1  34.00 0.58045  
husbandry:time 1.4337  1  34.00 0.23945  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.igfbp2)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.igfbp2,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp2.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response   SE   df lower.CL upper.CL
 conventional      663 73.5 11.4      520      846
 ecological        447 48.6 12.7      353      566

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio   SE df null t.ratio p.value
 conventional / ecological  1.48 0.23 12    1   2.540  0.0259

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp2),
         edf = df.residual(mod.igfbp2))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)     1.63 0.657 12    0.198     3.06

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.2422 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.igfbp2,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp2.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response   SE   df lower.CL upper.CL
 105dpc      536 44.3 15.4      450      639
 8dpp        553 45.8 15.4      464      660

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.969 0.0553 34    1  -0.558  0.5805

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp2),
         edf = df.residual(mod.igfbp2))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)   -0.132 0.236 34   -0.611    0.348

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.2422 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.igfbp2,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp2.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response   SE   df lower.CL upper.CL
 conventional      676 79.7 14.5      525      869
 ecological        425 49.3 16.4      333      544

time = 8dpp:
 husbandry    response   SE   df lower.CL upper.CL
 conventional      651 76.8 14.5      506      838
 ecological        470 54.5 16.4      368      601

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological  1.59 0.263 15.4    1   2.797  0.0133

time = 8dpp:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological  1.39 0.229 15.4    1   1.971  0.0670

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp2),
         edf = df.residual(mod.igfbp2))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)     1.91 0.703 15.4    0.415     3.41

time = 8dpp:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)     1.35 0.693 15.4   -0.128     2.82

sigma used for effect sizes: 0.2422 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.igfbp2,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp2.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response   SE   df lower.CL upper.CL
 105dpc      676 79.7 14.5      525      869
 8dpp        651 76.8 14.5      506      838

husbandry = ecological:
 time   response   SE   df lower.CL upper.CL
 105dpc      425 49.3 16.4      333      544
 8dpp        470 54.5 16.4      368      601

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 1.037 0.0837 34    1   0.452  0.6541

husbandry = ecological:
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.905 0.0730 34    1  -1.241  0.2230

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp2),
         edf = df.residual(mod.igfbp2))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    0.151 0.334 34   -0.527    0.829

husbandry = ecological:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)   -0.414 0.335 34   -1.095    0.268

sigma used for effect sizes: 0.2422 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.igfbp2 <- dat.l  |> 
  dplyr::filter(parameter == "igfbp2.ser") |> 
  droplevels() |> 
  drop_na(value) |> 
  # plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 1250), 
                     breaks = seq(0, 1250, 250),
                     minor_breaks = seq(0, 1250, by = 100) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "IGF BP2 (serum) [ng/ml]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.igfbp2

IGFBP3 (serum)

Model

hist(dat.l |> 
       dplyr::filter(parameter == "igfbp3.ser") |> 
       droplevels() |> 
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
         dplyr::filter(parameter == "igfbp3.ser") |> 
         droplevels() |> 
         drop_na(value) |> 
         pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.igfbp3 <- lmerTest::lmer(log(value) ~ 
                               husbandry +
                               time +
                               # interaction term
                               husbandry : time +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "igfbp3.ser") |> 
                               droplevels() |> 
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.igfbp3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp3.ser")),  
    value)
Control: contr

REML criterion at convergence: 159.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.78859 -0.50549  0.09757  0.46868  1.52002 

Random effects:
 Groups        Name        Variance  Std.Dev. 
 litter.no:sow (Intercept) 5.588e-01 7.475e-01
 sow           (Intercept) 6.984e-16 2.643e-08
 Residual                  2.038e-01 4.515e-01
Number of obs: 72, groups:  litter.no:sow, 36; sow, 15

Fixed effects:
                             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                   6.72944    0.20583 44.24476  32.694  < 2e-16 ***
husbandryecological           0.08475    0.29109 44.24476   0.291  0.77229    
time8dpp                      0.53588    0.15048 34.00000   3.561  0.00112 ** 
husbandryecological:time8dpp -0.03664    0.21282 34.00000  -0.172  0.86434    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.707              
time8dpp    -0.366  0.258       
hsbndrycl:8  0.258 -0.366 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.igfbp3, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time  0.006   0.006     1    34    0.03  0.864

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.igfbp3,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                 Chisq Df Pr(>Chisq)    
husbandry       0.0601  1     0.8063    
time           23.6578  1  1.151e-06 ***
husbandry:time  0.0296  1     0.8633    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.igfbp3,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                     F Df Df.res    Pr(>F)    
husbandry       0.0583  1 10.995    0.8136    
time           23.6578  1 34.000 2.577e-05 ***
husbandry:time  0.0296  1 34.000    0.8643    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.igfbp3)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.igfbp3,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp3.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response  SE df lower.CL upper.CL
 conventional     1094 210 34      741     1614
 ecological       1169 224 34      792     1725

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio    SE df null t.ratio p.value
 conventional / ecological 0.936 0.254 34    1  -0.245  0.8078

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp3),
         edf = df.residual(mod.igfbp3))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate  SE df lower.CL upper.CL
 (conventional - ecological)   -0.147 0.6 34    -1.37     1.07

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.4515 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.igfbp3,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp3.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response  SE   df lower.CL upper.CL
 105dpc      873 127 44.2      651     1170
 8dpp       1465 213 44.2     1092     1964

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.596 0.0634 34    1  -4.864  <.0001

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp3),
         edf = df.residual(mod.igfbp3))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    -1.15 0.256 34    -1.67   -0.626

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.4515 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.igfbp3,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp3.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response  SE   df lower.CL upper.CL
 conventional      837 172 44.2      553     1267
 ecological        911 187 44.2      602     1379

time = 8dpp:
 husbandry    response  SE   df lower.CL upper.CL
 conventional     1430 294 44.2      944     2165
 ecological       1500 309 44.2      991     2271

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological 0.919 0.267 44.2    1  -0.291  0.7723

time = 8dpp:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological 0.953 0.277 44.2    1  -0.165  0.8695

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp3),
         edf = df.residual(mod.igfbp3))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)   -0.188 0.645 44.2    -1.49     1.11

time = 8dpp:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)   -0.107 0.645 44.2    -1.41     1.19

sigma used for effect sizes: 0.4515 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.igfbp3,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp3.ser") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response  SE   df lower.CL upper.CL
 105dpc      837 172 44.2      553     1267
 8dpp       1430 294 44.2      944     2165

husbandry = ecological:
 time   response  SE   df lower.CL upper.CL
 105dpc      911 187 44.2      602     1379
 8dpp       1500 309 44.2      991     2271

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.585 0.0881 34    1  -3.561  0.0011

husbandry = ecological:
 contrast      ratio     SE df null t.ratio p.value
 105dpc / 8dpp 0.607 0.0913 34    1  -3.318  0.0022

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp3),
         edf = df.residual(mod.igfbp3))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    -1.19 0.349 34    -1.90   -0.477

husbandry = ecological:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    -1.11 0.347 34    -1.81   -0.400

sigma used for effect sizes: 0.4515 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.igfbp3 <- dat.l  |> 
  dplyr::filter(parameter == "igfbp3.ser") |> 
  droplevels() |> 
  drop_na(value) |> 
  # plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 6500), 
                     breaks = seq(0, 6500, 2000),
                     minor_breaks = seq(0, 6500, by = 200) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "IGF BP3 (serum) [ng/ml]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.igfbp3

IGFBP2/IGFBP3 (serum)

Model

hist(dat.l |> 
       dplyr::filter(parameter == "igfbp23.ser") |> 
       droplevels() |> 
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
         dplyr::filter(parameter == "igfbp23.ser") |> 
         droplevels() |>
         drop_na(value) |> 
         pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.igfbp23 <- lmerTest::lmer(log(value) ~ 
                               husbandry +
                               time +
                               # interaction term
                               husbandry : time +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "igfbp23.ser") |> 
                               droplevels() |>
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.igfbp23)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp23.ser")),  
    value)
Control: contr

REML criterion at convergence: 144.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.97572 -0.51849 -0.00712  0.49478  1.73489 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.4515   0.6719  
 sow           (Intercept) 0.0000   0.0000  
 Residual                  0.1608   0.4010  
Number of obs: 72, groups:  litter.no:sow, 36; sow, 15

Fixed effects:
                             Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)                   -0.2861     0.1844 44.0488  -1.551 0.127980    
husbandryecological           -0.5229     0.2608 44.0488  -2.005 0.051155 .  
time8dpp                      -0.5724     0.1337 34.0000  -4.282 0.000143 ***
husbandryecological:time8dpp   0.1733     0.1890 34.0000   0.917 0.365645    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.707              
time8dpp    -0.362  0.256       
hsbndrycl:8  0.256 -0.362 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.igfbp23, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time  0.135   0.135     1    34   0.841  0.366

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.igfbp23,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                 Chisq Df Pr(>Chisq)    
husbandry       3.2202  1    0.07273 .  
time           26.4099  1  2.761e-07 ***
husbandry:time  0.8407  1    0.35919    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.igfbp23,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                     F Df Df.res    Pr(>F)    
husbandry       3.1232  1 10.995    0.1049    
time           26.4099  1 34.000 1.134e-05 ***
husbandry:time  0.8407  1 34.000    0.3656    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.igfbp23)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.igfbp23,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp23.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response     SE df lower.CL upper.CL
 conventional    0.564 0.0970 34    0.398    0.800
 ecological      0.365 0.0627 34    0.257    0.517

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio    SE df null t.ratio p.value
 conventional / ecological  1.55 0.376 34    1   1.795  0.0816

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp23),
         edf = df.residual(mod.igfbp23))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)     1.09 0.614 34   -0.159     2.34

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.401 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.igfbp23,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp23.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response     SE df lower.CL upper.CL
 105dpc    0.578 0.0754 44    0.445    0.752
 8dpp      0.356 0.0464 44    0.274    0.463

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp  1.63 0.154 34    1   5.139  <.0001

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp23),
         edf = df.residual(mod.igfbp23))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)     1.21 0.259 34    0.686     1.74

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.401 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.igfbp23,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp23.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response     SE df lower.CL upper.CL
 conventional    0.751 0.1390 44    0.518    1.089
 ecological      0.445 0.0821 44    0.307    0.646

time = 8dpp:
 husbandry    response     SE df lower.CL upper.CL
 conventional    0.424 0.0782 44    0.292    0.615
 ecological      0.299 0.0551 44    0.206    0.433

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio   SE df null t.ratio p.value
 conventional / ecological  1.69 0.44 44    1   2.005  0.0512

time = 8dpp:
 contrast                  ratio   SE df null t.ratio p.value
 conventional / ecological  1.42 0.37 44    1   1.340  0.1870

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp23),
         edf = df.residual(mod.igfbp23))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)    1.304 0.660 44  -0.0269     2.64

time = 8dpp:
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)    0.872 0.655 44  -0.4481     2.19

sigma used for effect sizes: 0.401 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.igfbp23,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "igfbp23.ser") |> 
          droplevels() |>
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response     SE df lower.CL upper.CL
 105dpc    0.751 0.1390 44    0.518    1.089
 8dpp      0.424 0.0782 44    0.292    0.615

husbandry = ecological:
 time   response     SE df lower.CL upper.CL
 105dpc    0.445 0.0821 44    0.307    0.646
 8dpp      0.299 0.0551 44    0.206    0.433

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp  1.77 0.237 34    1   4.282  0.0001

husbandry = ecological:
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp  1.49 0.199 34    1   2.986  0.0052

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.igfbp23),
         edf = df.residual(mod.igfbp23))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    1.427 0.356 34    0.704     2.15

husbandry = ecological:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    0.995 0.345 34    0.295     1.70

sigma used for effect sizes: 0.401 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.igfbp23 <- dat.l  |> 
  dplyr::filter(parameter == "igfbp23.ser") |> 
  droplevels() |>
  drop_na(value) |> 
  # plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = (value))) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 2.5), 
                     breaks = seq(0, 2.5, 0.5),
                     minor_breaks = seq(0, 2.5, by = 0.1) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "Molar ratio of IGFBP2/IGFBP3 (serum)",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.igfbp23

Proteolytic activity (serum)

Model

hist(dat.l |> 
       dplyr::filter(parameter == "proteolysis") |> 
       droplevels() |> 
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
         dplyr::filter(parameter == "proteolysis") |> 
         droplevels() |> 
         drop_na(value) |> 
         pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.prot <- lmerTest::lmer(log(value+1) ~ 
                               husbandry +
                               time +
                               # interaction term
                               husbandry : time +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                              dplyr::filter(parameter == "proteolysis") |> 
                              droplevels() |> 
                              drop_na(value),
                             REML    = TRUE,
                             control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.prot)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
log(value + 1) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "proteolysis")),  
    value)
Control: contr

REML criterion at convergence: 120.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4098 -1.1650  0.3079  0.7632  1.3081 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.000    0.000   
 sow           (Intercept) 0.000    0.000   
 Residual                  1.413    1.189   
Number of obs: 39, groups:  litter.no:sow, 20; sow, 15

Fixed effects:
                             Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)                    1.6496     0.3963 35.0000   4.163 0.000194 ***
husbandryecological           -0.2007     0.5462 35.0000  -0.367 0.715542    
time8dpp                      -0.2647     0.5462 35.0000  -0.485 0.630926    
husbandryecological:time8dpp   0.4917     0.7622 35.0000   0.645 0.523064    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.725              
time8dpp    -0.725  0.526       
hsbndrycl:8  0.520 -0.717 -0.717
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.prot, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value + 1) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time  0.588   0.588     1    35   0.416  0.523

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.prot,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value + 1)
                Chisq Df Pr(>Chisq)
husbandry      0.0185  1     0.8918
time           0.0010  1     0.9744
husbandry:time 0.4162  1     0.5189
car::Anova(mod.prot,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value + 1)
                    F Df Df.res Pr(>F)
husbandry      0.0175  1 10.932 0.8972
time           0.0007  1 17.691 0.9796
husbandry:time 0.4121  1 17.726 0.5291

Model diagnostics

performance::check_model(mod.prot)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.prot,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "proteolysis") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response   SE df lower.CL upper.CL
 conventional     3.56 1.25 35     1.62     6.94
 ecological       3.77 1.27 35     1.78     7.18

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log(mu + 1) scale 

$contrasts
 contrast                  estimate    SE df t.ratio p.value
 conventional - ecological  -0.0452 0.381 35  -0.119  0.9063

Results are averaged over the levels of: time 
Note: contrasts are still on the log(mu + 1) scale. Consider using
      regrid() if you want contrasts of back-transformed estimates. 
Degrees-of-freedom method: satterthwaite 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.prot),
         edf = df.residual(mod.prot))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)   -0.038 0.321 35   -0.689    0.613

Results are averaged over the levels of: time 
sigma used for effect sizes: 1.189 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.prot,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "proteolysis") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response   SE df lower.CL upper.CL
 105dpc     3.71 1.29 35     1.70     7.20
 8dpp       3.62 1.23 35     1.69     6.93

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log(mu + 1) scale 

$contrasts
 contrast      estimate    SE df t.ratio p.value
 105dpc - 8dpp   0.0189 0.381 35   0.050  0.9608

Results are averaged over the levels of: husbandry 
Note: contrasts are still on the log(mu + 1) scale. Consider using
      regrid() if you want contrasts of back-transformed estimates. 
Degrees-of-freedom method: satterthwaite 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.prot),
         edf = df.residual(mod.prot))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)   0.0159 0.321 35   -0.635    0.667

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 1.189 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.prot,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "proteolysis") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response   SE df lower.CL upper.CL
 conventional     4.21 2.06 35    1.328    10.64
 ecological       3.26 1.60 35    0.985     8.14

time = 8dpp:
 husbandry    response   SE df lower.CL upper.CL
 conventional     2.99 1.50 35    0.862     7.57
 ecological       4.34 2.01 35    1.491    10.46

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log(mu + 1) scale 

$contrasts
time = 105dpc:
 contrast                  estimate    SE df t.ratio p.value
 conventional - ecological    0.201 0.546 35   0.367  0.7155

time = 8dpp:
 contrast                  estimate    SE df t.ratio p.value
 conventional - ecological   -0.291 0.532 35  -0.547  0.5875

Note: contrasts are still on the log(mu + 1) scale. Consider using
      regrid() if you want contrasts of back-transformed estimates. 
Degrees-of-freedom method: satterthwaite 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.prot),
         edf = df.residual(mod.prot))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)    0.169 0.460 35   -0.765    1.103

time = 8dpp:
 contrast                    estimate    SE df lower.CL upper.CL
 (conventional - ecological)   -0.245 0.448 35   -1.155    0.665

sigma used for effect sizes: 1.189 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.prot,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "proteolysis") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response   SE df lower.CL upper.CL
 105dpc     4.21 2.06 35    1.328    10.64
 8dpp       2.99 1.50 35    0.862     7.57

husbandry = ecological:
 time   response   SE df lower.CL upper.CL
 105dpc     3.26 1.60 35    0.985     8.14
 8dpp       4.34 2.01 35    1.491    10.46

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log(mu + 1) scale 

$contrasts
husbandry = conventional:
 contrast      estimate    SE df t.ratio p.value
 105dpc - 8dpp    0.265 0.546 35   0.485  0.6309

husbandry = ecological:
 contrast      estimate    SE df t.ratio p.value
 105dpc - 8dpp   -0.227 0.532 35  -0.427  0.6720

Note: contrasts are still on the log(mu + 1) scale. Consider using
      regrid() if you want contrasts of back-transformed estimates. 
Degrees-of-freedom method: satterthwaite 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.prot),
         edf = df.residual(mod.prot))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    0.223 0.460 35   -0.712    1.157

husbandry = ecological:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)   -0.191 0.448 35   -1.100    0.718

sigma used for effect sizes: 1.189 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.prot <- dat.l  |> 
  dplyr::filter(parameter == "proteolysis") |> 
  droplevels() |> 
  drop_na(value) |> 
  # plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 21), 
                     breaks = seq(0, 21, 5),
                     minor_breaks = seq(0, 21, by = 1) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "Proteolytic activity (serum) [%]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.prot

Stanniocalcin, STC1 (serum)

Model

hist(dat.l |> 
       dplyr::filter(parameter == "stc1") |> 
       droplevels() |> 
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
         dplyr::filter(parameter == "stc1") |> 
         droplevels() |> 
         drop_na(value) |> 
         pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.stc1 <- lmerTest::lmer(log(value) ~ 
                               husbandry +
                               time +
                               # interaction term
                               husbandry : time +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "stc1") |> 
                               droplevels() |> 
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.stc1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "stc1")),  
    value)
Control: contr

REML criterion at convergence: 88.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.58308 -0.48078  0.00956  0.47032  1.39341 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 3.17571  1.7821  
 sow           (Intercept) 0.00000  0.0000  
 Residual                  0.05628  0.2372  
Number of obs: 38, groups:  litter.no:sow, 19; sow, 14

Fixed effects:
                             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                   5.62613    0.59926 17.29863   9.388 3.32e-08 ***
husbandryecological           1.99293    0.82602 17.29863   2.413   0.0272 *  
time8dpp                     -0.31643    0.11184 17.00000  -2.829   0.0116 *  
husbandryecological:time8dpp  0.08471    0.15416 17.00000   0.550   0.5898    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.725              
time8dpp    -0.093  0.068       
hsbndrycl:8  0.068 -0.093 -0.725
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.stc1, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time  0.017   0.017     1    17   0.302   0.59

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.stc1,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                 Chisq Df Pr(>Chisq)    
husbandry       6.1244  1  0.0133326 *  
time           12.4733  1  0.0004128 ***
husbandry:time  0.3020  1  0.5826404    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.stc1,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                     F Df Df.res   Pr(>F)   
husbandry       5.4881  1 10.123 0.040863 * 
time           12.4733  1 17.000 0.002561 **
husbandry:time  0.3020  1 17.000 0.589787   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.stc1)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.stc1,
        pairwise ~ husbandry, 
        data    = dat.l |>
          dplyr::filter(parameter == "stc1") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response   SE df lower.CL upper.CL
 conventional      237  141 17     67.3      834
 ecological       1814 1030 17    549.5     5987

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio    SE df null t.ratio p.value
 conventional / ecological 0.131 0.107 17    1  -2.475  0.0242

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.stc1),
         edf = df.residual(mod.stc1))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate   SE df lower.CL upper.CL
 (conventional - ecological)    -8.58 3.63 17    -16.2   -0.912

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.2372 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.stc1,
        pairwise ~ time, 
        data    = dat.l |>
          dplyr::filter(parameter == "stc1") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response  SE   df lower.CL upper.CL
 105dpc      752 311 17.3      315     1795
 8dpp        572 236 17.3      239     1365

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp  1.32 0.101 17    1   3.556  0.0024

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.stc1),
         edf = df.residual(mod.stc1))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)     1.16 0.356 17    0.403     1.91

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.2372 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.stc1,
        pairwise ~ husbandry|time, 
        data    = dat.l |>
          dplyr::filter(parameter == "stc1") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response   SE   df lower.CL upper.CL
 conventional      278  166 17.3     78.5      981
 ecological       2037 1160 17.3    614.7     6747

time = 8dpp:
 husbandry    response   SE   df lower.CL upper.CL
 conventional      202  121 17.3     57.2      715
 ecological       1615  918 17.3    487.6     5352

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological 0.136 0.113 17.3    1  -2.413  0.0272

time = 8dpp:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological 0.125 0.103 17.3    1  -2.515  0.0220

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.stc1),
         edf = df.residual(mod.stc1))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate   SE   df lower.CL upper.CL
 (conventional - ecological)    -8.40 3.64 17.3    -16.1   -0.727

time = 8dpp:
 contrast                    estimate   SE   df lower.CL upper.CL
 (conventional - ecological)    -8.76 3.66 17.3    -16.5   -1.056

sigma used for effect sizes: 0.2372 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.stc1,
        pairwise ~ time|husbandry, 
        data    = dat.l |>
          dplyr::filter(parameter == "stc1") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response   SE   df lower.CL upper.CL
 105dpc      278  166 17.3     78.5      981
 8dpp        202  121 17.3     57.2      715

husbandry = ecological:
 time   response   SE   df lower.CL upper.CL
 105dpc     2037 1160 17.3    614.7     6747
 8dpp       1615  918 17.3    487.6     5352

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp  1.37 0.153 17    1   2.829  0.0116

husbandry = ecological:
 contrast      ratio    SE df null t.ratio p.value
 105dpc / 8dpp  1.26 0.134 17    1   2.184  0.0433

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.stc1),
         edf = df.residual(mod.stc1))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    1.334 0.501 17  0.27695     2.39

husbandry = ecological:
 contrast        estimate    SE df lower.CL upper.CL
 (105dpc - 8dpp)    0.977 0.464 17 -0.00245     1.96

sigma used for effect sizes: 0.2372 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.stc1 <- dat.l  |> 
  dplyr::filter(parameter == "stc1") |> 
  droplevels() |> 
  drop_na(value) |> 
  # plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 10000), 
                     breaks = seq(0, 10000, 2000),
                     minor_breaks = seq(0, 10000, by = 400) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "STC1 (serum) [pg/ml]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.stc1

Combined plots: Figure 6

# Combine plots with a designated area for the legend
combined <- (plot.igfbp23 + 
             plot.stc1) + 
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(face = "bold", size = 20),
        legend.position = "bottom")
combined

png("./plots/figure6.png",
     width = 300, height = 200, units = "mm",
     pointsize = 10, res = 600)

combined

dev.off()
png 
  2 

Calcium (serum)

Model

hist(dat.l |> 
       dplyr::filter(parameter == "calc") |> 
       droplevels() |> 
       drop_na(value) |> 
       pull("value"), 
     breaks = 30)

hist(log(dat.l |> 
         dplyr::filter(parameter == "calc") |> 
         droplevels() |> 
         drop_na(value) |> 
         pull("value")),
     breaks = 30)

Singularity issues if nested structure defined like (1 | sow/litter.no).

mod.calc <- lmerTest::lmer(log(value) ~ 
                               husbandry +
                               time +
                               # interaction term
                               husbandry : time +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "calc") |> 
                               dplyr::filter(time != "30dpc") |> 
                               droplevels() |> 
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.calc)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dplyr::filter(dat.l, parameter ==  
    "calc"), time != "30dpc")), value)
Control: contr

REML criterion at convergence: 4.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.9530 -0.3184  0.0786  0.3944  1.6808 

Random effects:
 Groups        Name        Variance  Std.Dev. 
 litter.no:sow (Intercept) 6.379e-18 2.526e-09
 sow           (Intercept) 3.982e-02 1.995e-01
 Residual                  3.860e-02 1.965e-01
Number of obs: 80, groups:  litter.no:sow, 40; sow, 15

Fixed effects:
                             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                   2.42516    0.08737 15.46761  27.756 1.29e-14 ***
husbandryecological          -0.08834    0.12130 16.18553  -0.728    0.477    
time8dpp                     -0.02237    0.06213 62.42643  -0.360    0.720    
husbandryecological:time8dpp  0.10889    0.08786 62.42643   1.239    0.220    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) hsbndr tm8dpp
hsbndryclgc -0.720              
time8dpp    -0.356  0.256       
hsbndrycl:8  0.251 -0.362 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.calc, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
               Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
husbandry:time  0.059   0.059     1 62.426   1.536   0.22

Without significant interactions, choose type = 2. With significant interactions, choose type = 3.

car::Anova(mod.calc,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
                Chisq Df Pr(>Chisq)
husbandry      0.0899  1     0.7644
time           0.5333  1     0.4652
husbandry:time 1.5360  1     0.2152
car::Anova(mod.calc,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
                    F Df Df.res Pr(>F)
husbandry      0.0897  1 12.849 0.7694
time           0.5333  1 38.000 0.4697
husbandry:time 1.5360  1 38.000 0.2228

Model diagnostics

performance::check_model(mod.calc)

Emmeans & Effect sizes

Husbandry

Emmeans:

emm <- emmeans(mod.calc,
        pairwise ~ husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "calc") |> 
          dplyr::filter(time != "30dpc") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 husbandry    response    SE   df lower.CL upper.CL
 conventional     11.2 0.913 11.8     9.35     13.4
 ecological       10.8 0.845 12.7     9.12     12.8

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological  1.03 0.117 12.2    1   0.300  0.7694

Results are averaged over the levels of: time 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.calc),
         edf = df.residual(mod.calc))
Since 'object' is a list, we are using the contrasts already present.
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)    0.173 0.576 12.2    -1.08     1.42

Results are averaged over the levels of: time 
sigma used for effect sizes: 0.1965 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point

Emmeans:

emm <- emmeans(mod.calc,
        pairwise ~ time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "calc") |> 
          dplyr::filter(time != "30dpc") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")
NOTE: Results may be misleading due to involvement in interactions
emm
$emmeans
 time   response    SE   df lower.CL upper.CL
 105dpc     10.8 0.656 16.2     9.51     12.3
 8dpp       11.2 0.677 16.2     9.82     12.7

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast      ratio     SE   df null t.ratio p.value
 105dpc / 8dpp 0.968 0.0425 62.4    1  -0.730  0.4680

Results are averaged over the levels of: husbandry 
Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.calc),
         edf = df.residual(mod.calc))
Since 'object' is a list, we are using the contrasts already present.
 contrast        estimate    SE   df lower.CL upper.CL
 (105dpc - 8dpp)   -0.163 0.224 62.4   -0.611    0.284

Results are averaged over the levels of: husbandry 
sigma used for effect sizes: 0.1965 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Husbandry:Time point

Emmeans:

emm <- emmeans(mod.calc,
        pairwise ~ husbandry|time, 
        data    = dat.l |> 
          dplyr::filter(parameter == "calc") |> 
          dplyr::filter(time != "30dpc") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
time = 105dpc:
 husbandry    response    SE   df lower.CL upper.CL
 conventional     11.3 0.988 15.5     9.39     13.6
 ecological       10.3 0.871 17.0     8.67     12.4

time = 8dpp:
 husbandry    response    SE   df lower.CL upper.CL
 conventional     11.1 0.966 15.5     9.18     13.3
 ecological       11.3 0.949 17.0     9.45     13.5

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
time = 105dpc:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological  1.09 0.133 16.2    1   0.728  0.4768

time = 8dpp:
 contrast                  ratio    SE   df null t.ratio p.value
 conventional / ecological  0.98 0.119 16.2    1  -0.169  0.8675

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.calc),
         edf = df.residual(mod.calc))
Since 'object' is a list, we are using the contrasts already present.
time = 105dpc:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)    0.450 0.619 16.2    -0.86     1.76

time = 8dpp:
 contrast                    estimate    SE   df lower.CL upper.CL
 (conventional - ecological)   -0.105 0.617 16.2    -1.41     1.20

sigma used for effect sizes: 0.1965 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Time point:Husbandry

Emmeans:

emm <- emmeans(mod.calc,
        pairwise ~ time|husbandry, 
        data    = dat.l |> 
          dplyr::filter(parameter == "calc") |> 
          dplyr::filter(time != "30dpc") |> 
          droplevels() |> 
          drop_na(value), 
        adjust  = "tukey",
        lmer.df = "satterthwaite",
        type    = "response")

emm
$emmeans
husbandry = conventional:
 time   response    SE   df lower.CL upper.CL
 105dpc     11.3 0.988 15.5     9.39     13.6
 8dpp       11.1 0.966 15.5     9.18     13.3

husbandry = ecological:
 time   response    SE   df lower.CL upper.CL
 105dpc     10.3 0.871 17.0     8.67     12.4
 8dpp       11.3 0.949 17.0     9.45     13.5

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
husbandry = conventional:
 contrast      ratio     SE   df null t.ratio p.value
 105dpc / 8dpp 1.023 0.0635 62.4    1   0.360  0.7201

husbandry = ecological:
 contrast      ratio     SE   df null t.ratio p.value
 105dpc / 8dpp 0.917 0.0570 62.4    1  -1.393  0.1686

Degrees-of-freedom method: satterthwaite 
Tests are performed on the log scale 

Effect sizes:

eff_size(emm,
         sigma = sigma(mod.calc),
         edf = df.residual(mod.calc))
Since 'object' is a list, we are using the contrasts already present.
husbandry = conventional:
 contrast        estimate    SE   df lower.CL upper.CL
 (105dpc - 8dpp)    0.114 0.316 62.4   -0.518    0.746

husbandry = ecological:
 contrast        estimate    SE   df lower.CL upper.CL
 (105dpc - 8dpp)   -0.440 0.318 62.4   -1.077    0.196

sigma used for effect sizes: 0.1965 
Degrees-of-freedom method: inherited from satterthwaite when re-gridding 
Confidence level used: 0.95 

Plot

plot.calc <- dat.l  |> 
  dplyr::filter(parameter == "calc") |> 
  dplyr::filter(time != "30dpc") |> 
  droplevels() |> 
  drop_na(value) |> 
  # plot
  mutate(jit = jitter(as.numeric(time), 0.3)) |>  
  ggplot(aes(y   = value)) +
  geom_boxplot(aes(x   = time, 
                   fill = husbandry),
               col = "black",
               outlier.shape = NA, width = 0.5) +
  geom_jitter(aes(x   = time,
                  fill = husbandry),
              col = "black",
              position = position_jitterdodge(jitter.width  = 0.15,
                                              dodge.width   = 0.5),
              size = 2, alpha = 0.7,
              show.legend = FALSE) +
  scale_fill_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 21), 
                     breaks = seq(0, 20, 5),
                     minor_breaks = seq(0, 20, by = 2) ) +
  scale_x_discrete(labels= c("105dpc", "8dpp")) +
  labs(x = "",
       y = "Calcium (serum) [mg/dl]",
       fill = "Husbandry") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         fill = guide_legend(override.aes = list(linetype = 0, size  = 5)))
plot.calc

How to cite R

“All analyses were performed using R Statistical Software (version 4.4.2; R Core Team 2024)”.

Reference: R Core Team (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

citation()
To cite R in publications use:

  R Core Team (2024). _R: A Language and Environment for Statistical
  Computing_. R Foundation for Statistical Computing, Vienna, Austria.
  <https://www.R-project.org/>.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Manual{,
    title = {R: A Language and Environment for Statistical Computing},
    author = {{R Core Team}},
    organization = {R Foundation for Statistical Computing},
    address = {Vienna, Austria},
    year = {2024},
    url = {https://www.R-project.org/},
  }

We have invested a lot of time and effort in creating R, please cite it
when using it for data analysis. See also 'citation("pkgname")' for
citing R packages.
version$version.string
[1] "R version 4.4.2 (2024-10-31 ucrt)"
citation("tidyverse")
Um Paket 'tidyverse' in Publikationen zu zitieren, nutzen Sie bitte:

  Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R,
  Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller
  E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V,
  Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to
  the tidyverse." _Journal of Open Source Software_, *4*(43), 1686.
  doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Article{,
    title = {Welcome to the {tidyverse}},
    author = {Hadley Wickham and Mara Averick and Jennifer Bryan and Winston Chang and Lucy D'Agostino McGowan and Romain François and Garrett Grolemund and Alex Hayes and Lionel Henry and Jim Hester and Max Kuhn and Thomas Lin Pedersen and Evan Miller and Stephan Milton Bache and Kirill Müller and Jeroen Ooms and David Robinson and Dana Paige Seidel and Vitalie Spinu and Kohske Takahashi and Davis Vaughan and Claus Wilke and Kara Woo and Hiroaki Yutani},
    year = {2019},
    journal = {Journal of Open Source Software},
    volume = {4},
    number = {43},
    pages = {1686},
    doi = {10.21105/joss.01686},
  }
citation("kableExtra")
Um Paket 'kableExtra' in Publikationen zu zitieren, nutzen Sie bitte:

  Zhu H (2024). _kableExtra: Construct Complex Table with 'kable' and
  Pipe Syntax_. R package version 1.4.0,
  <https://CRAN.R-project.org/package=kableExtra>.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Manual{,
    title = {kableExtra: Construct Complex Table with 'kable' and Pipe Syntax},
    author = {Hao Zhu},
    year = {2024},
    note = {R package version 1.4.0},
    url = {https://CRAN.R-project.org/package=kableExtra},
  }
citation("patchwork")
Um Paket 'patchwork' in Publikationen zu zitieren, nutzen Sie bitte:

  Pedersen T (2024). _patchwork: The Composer of Plots_. R package
  version 1.3.0, <https://CRAN.R-project.org/package=patchwork>.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Manual{,
    title = {patchwork: The Composer of Plots},
    author = {Thomas Lin Pedersen},
    year = {2024},
    note = {R package version 1.3.0},
    url = {https://CRAN.R-project.org/package=patchwork},
  }
citation("lmerTest")
To cite lmerTest in publications use:

  Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest
  Package: Tests in Linear Mixed Effects Models." _Journal of
  Statistical Software_, *82*(13), 1-26. doi:10.18637/jss.v082.i13
  <https://doi.org/10.18637/jss.v082.i13>.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Article{,
    title = {{lmerTest} Package: Tests in Linear Mixed Effects Models},
    author = {Alexandra Kuznetsova and Per B. Brockhoff and Rune H. B. Christensen},
    journal = {Journal of Statistical Software},
    year = {2017},
    volume = {82},
    number = {13},
    pages = {1--26},
    doi = {10.18637/jss.v082.i13},
  }
citation("car")
To cite the car package in publications use:

  Fox J, Weisberg S (2019). _An R Companion to Applied Regression_,
  Third edition. Sage, Thousand Oaks CA.
  <https://www.john-fox.ca/Companion/>.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Book{,
    title = {An {R} Companion to Applied Regression},
    edition = {Third},
    author = {John Fox and Sanford Weisberg},
    year = {2019},
    publisher = {Sage},
    address = {Thousand Oaks {CA}},
    url = {https://www.john-fox.ca/Companion/},
  }
citation("emmeans")
Um Paket 'emmeans' in Publikationen zu zitieren, nutzen Sie bitte:

  Lenth R (2024). _emmeans: Estimated Marginal Means, aka Least-Squares
  Means_. R package version 1.10.6,
  <https://CRAN.R-project.org/package=emmeans>.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Manual{,
    title = {emmeans: Estimated Marginal Means, aka Least-Squares Means},
    author = {Russell V. Lenth},
    year = {2024},
    note = {R package version 1.10.6},
    url = {https://CRAN.R-project.org/package=emmeans},
  }
citation("performance")
Um Paket 'performance' in Publikationen zu zitieren, nutzen Sie bitte:

  Lüdecke et al., (2021). performance: An R Package for Assessment,
  Comparison and Testing of Statistical Models. Journal of Open Source
  Software, 6(60), 3139. https://doi.org/10.21105/joss.03139

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Article{,
    title = {{performance}: An {R} Package for Assessment, Comparison and Testing of Statistical Models},
    author = {Daniel Lüdecke and Mattan S. Ben-Shachar and Indrajeet Patil and Philip Waggoner and Dominique Makowski},
    year = {2021},
    journal = {Journal of Open Source Software},
    volume = {6},
    number = {60},
    pages = {3139},
    doi = {10.21105/joss.03139},
  }
citation("Hmisc")
Um Paket 'Hmisc' in Publikationen zu zitieren, nutzen Sie bitte:

  Harrell Jr F (2025). _Hmisc: Harrell Miscellaneous_. R package
  version 5.2-2, <https://CRAN.R-project.org/package=Hmisc>.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Manual{,
    title = {Hmisc: Harrell Miscellaneous},
    author = {Frank E {Harrell Jr}},
    year = {2025},
    note = {R package version 5.2-2},
    url = {https://CRAN.R-project.org/package=Hmisc},
  }
citation("viridis")
To cite viridis/viridisLite in publications use:

  Simon Garnier, Noam Ross, Robert Rudis, Antônio P. Camargo, Marco
  Sciaini, and Cédric Scherer (2024). viridis(Lite) -
  Colorblind-Friendly Color Maps for R. viridis package version 0.6.5.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Manual{,
    title = {{viridis(Lite)} - Colorblind-Friendly Color Maps for R},
    author = {{Garnier} and {Simon} and {Ross} and {Noam} and {Rudis} and {Robert} and {Camargo} and Antônio Pedro and {Sciaini} and {Marco} and {Scherer} and {Cédric}},
    year = {2024},
    note = {viridis package version 0.6.5},
    url = {https://sjmgarnier.github.io/viridis/},
    doi = {10.5281/zenodo.4679423},
  }

Session Info

sessionInfo()
R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default


locale:
[1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8   
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] viridis_0.6.5      viridisLite_0.4.2  Hmisc_5.2-2        performance_0.13.0
 [5] emmeans_1.10.6     car_3.1-3          carData_3.0-5      lmerTest_3.1-3    
 [9] lme4_1.1-36        Matrix_1.7-1       patchwork_1.3.0    kableExtra_1.4.0  
[13] lubridate_1.9.4    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
[17] purrr_1.0.2        readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
[21] ggplot2_3.5.1      tidyverse_2.0.0   

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1    farver_2.1.2        fastmap_1.2.0      
 [4] TH.data_1.1-3       bayestestR_0.15.1   rpart_4.1.23       
 [7] digest_0.6.37       timechange_0.3.0    estimability_1.5.1 
[10] lifecycle_1.0.4     cluster_2.1.6       survival_3.7-0     
[13] magrittr_2.0.3      compiler_4.4.2      rlang_1.1.4        
[16] tools_4.4.2         yaml_2.3.10         data.table_1.16.4  
[19] knitr_1.49          labeling_0.4.3      htmlwidgets_1.6.4  
[22] xml2_1.3.6          abind_1.4-8         multcomp_1.4-26    
[25] withr_3.0.2         foreign_0.8-87      numDeriv_2016.8-1.1
[28] datawizard_1.0.0    nnet_7.3-19         grid_4.4.2         
[31] xtable_1.8-4        colorspace_2.1-1    scales_1.3.0       
[34] MASS_7.3-61         insight_1.0.1       cli_3.6.3          
[37] mvtnorm_1.3-3       rmarkdown_2.29      reformulas_0.4.0   
[40] generics_0.1.3      rstudioapi_0.17.1   tzdb_0.4.0         
[43] minqa_1.2.8         splines_4.4.2       parallel_4.4.2     
[46] base64enc_0.1-3     vctrs_0.6.5         boot_1.3-31        
[49] sandwich_3.1-1      jsonlite_1.8.9      hms_1.1.3          
[52] ggrepel_0.9.6       pbkrtest_0.5.3      htmlTable_2.4.3    
[55] Formula_1.2-5       systemfonts_1.1.0   see_0.9.0          
[58] glue_1.8.0          nloptr_2.1.1        codetools_0.2-20   
[61] stringi_1.8.4       gtable_0.3.6        munsell_0.5.1      
[64] pillar_1.10.1       htmltools_0.5.8.1   R6_2.6.1           
[67] Rdpack_2.6.2        evaluate_1.0.3      lattice_0.22-6     
[70] backports_1.5.0     rbibutils_2.3       broom_1.0.7        
[73] Rcpp_1.0.13-1       checkmate_2.3.2     gridExtra_2.3      
[76] svglite_2.1.3       coda_0.19-4.1       nlme_3.1-166       
[79] mgcv_1.9-1          xfun_0.50           zoo_1.8-12         
[82] pkgconfig_2.0.3